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The protein p53 is a key regulator of cellular response to a wide variety of stressors. In 
cancer cells inhibitory regulators of p53 such as MDM2 and MDMX proteins are often 
overexpressed. We apply in silico techniques to better understand the role and interac- 
tions of these proteins in a cell cycle process. Furthermore we investigate the role of 
stochasticity in determining system behavior. We have found that stochasticity is able to 
affect system behavior profoundly. We also derive a general result for the way in which ini- 
tially synchronized oscillating stochastic systems will fall out of synchronization with each 
other. 
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INTRODUCTION 

Among the vast number of mechanisms utilized by cancer cells to 
sustain cell division, the inactivation of the essential tumor sup- 
pressor and transcription factor p53 is one of the most frequent 
and effective strategies. Therefore, restoring the activity of the 
p53 -signaling pathway is currently one of the most promising ther- 
apeutic strategies for fighting this disease (Levine and Oren, 2009). 

In normal cells, p53 plays a central role in the regulation of the 
cell cycle, apoptosis, DNA repair, and senescence (Teodoro et al., 
2007); p53 responds to cellular stress, such as hypoxia or DNA 
damage, by accumulating in the nucleus, regulating the expression 
of target genes, and activating/ inactivating various pathways in 
order to maintain the normal function of the cell (Maltzman and 
Czyzyk, 1984; Kastan et al, 1991; Graeber et al, 1994). Indeed, 
it appears that whenever the integrity of a cell's genetic code is 
threatened, p53 is there to protect it. This conclusion has led p53 
to be called the guardian of the genome (Lane, 1992). 

However, the p53 -signaling pathway is inoperative in almost 
all types of human cancer; factors that inactivate p53 specifically 
include genetic mutations or deletions (Feki and Irminger- Finger, 
2004), defective post-translational modifications, and interactions 
with its main endogenous inhibitors, MDM2 (Momand et al., 
1998) and MDMX (Shvarts et al, 1996). Excitingly, a number 
of these tumors have been shown to have a less invasive phenotype 
upon restoration of p53 activity (Olivier et al., 2002; Ventura et al., 
2007; Suad et al, 2009; Mandinova and Lee, 201 1). 

With the cost of drug development on the scale of hundreds 
of millions to billions of dollars per new drug entity - and ris- 
ing - there is strong need to look for any possible acceleration and 
improvement to the efficiency and accuracy of the development 
process (Paul et al., 2010). Thanks to the increasing comput- 
ing power available to researchers, it is now becoming affordable 
and practical to attempt to use in silico models to improve the 



development process. One way to do this is to improve the ability of 
researchers to select appropriate proteins, or interactions between 
proteins, as targets for drug development by better understanding 
their function in protein interaction networks. 

The purpose of this study is to gain new insights into the func- 
tioning of p53, a central protein in cell cycle regulation. A simple 
model of p53 oscillations in response to ionizing radiation is pre- 
sented. Additionally, the behavior of stochastic and deterministic 
representations of the same model system is compared. 

CELL CYCLE 

The protein p53 is a regulator of the cell cycle and cell fate. Under 
normal conditions, a cell will normally progress through several 
stages. In the Gl phase (first gap phase) the cell grows in size to 
prepare for DNA synthesis. After Gl, the cell moves into S phase 
(synthesis phase), during which new DNA is synthesized. Cells that 
are not replicating can also leave Gl and enter the GO phase, a state 
in which they do not grow, and can remain quiescent indefinitely. 
Next comes the G2 phase (second gap phase), where cells grow 
further and complete their final preparations for mitosis. Mitosis 
then occurs and the cycle can begin anew (Lodish et al., 2008). A 
damaged cell may need to halt its cycle or even self-destruct in a 
process called apoptosis. Apoptosis is necessary for normal devel- 
opment and homeostasis of multicellular organisms, and is also a 
desirable outcome for cancer cells during cancer chemotherapy. 

In order to ensure that the process of cell division is carefully 
regulated, the cell has a number of checkpoints. These checkpoints 
are conditions that a cell must meet in order to progress in the 
cell cycle. For example, one checkpoint in Gl ensures that a cell 
has grown sufficiently in size to move into S phase and replicate 
its DNA. Another checkpoint that occurs in Gl is mediated by 
the protein p53: when DNA is damaged, p53 halts the cell cycle 
until the damage is repaired; this prevents the cell from trying 
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to duplicate the damaged DNA. When p53 is inactivated, this 
checkpoint no longer functions. A cell attempting to duplicate 
damaged DNA is likely to accumulate mutations (Alberts et al, 
1994). Figure 1 diagrams the relevance of p53 to the cell cycle. 

p53 

The protein p53 responds to many stressors including ultravio- 
let light (Maltzman and Czyzyk, 1984), ionizing radiation (Kastan 
et al., 1991), hypoxia, heat (Graeber et al., 1994), improper cell 
adhesion (Nigro et al., 1997), ribonucleotide depletion (Linke 
et al., 1996), and infection by influenza (Turpin et al., 2005). Some 
viral proteins are known to interact with p53, for example hepatitis 
B virus HBx protein (Truant et al., 1995) and the large T antigen 
of simian virus 40 (Dobbelstein and Roth, 1998). The protein p53 
has been demonstrated to induce cell cycle arrest, senescence, and 
apoptosis, with the specific outcome dependent on the extent and 
type of stress, and the genetic background of the cell (Vousden 
and Lu, 2002). The expression of p53 is tightly regulated by the 
cell (Sugrue et al., 1997; Lodish et al., 2008). In order to help it exe- 
cute its various functions p53 is post-translationally modified at 
many sites to determine its response (Meek and Anderson, 2009; 
Dai and Gu, 2010). The protein p53 transcriptionally regulates 
numerous genes, with a pattern that varies depending on the type 
of stress and the cell type (Zhao et al., 2000). In addition to its tran- 
scriptional activity, p53 plays a transcription-independent role in 
apoptosis by binding to several anti-apoptotic proteins (Mihara 
etal.,2003). 

The protein p53 is known to be mutated in approximately 50% 
of human tumors (Soussi and Wiman, 2007; Brown et al., 2009; 
Freed-Pastor and Prives, 2012). In addition, in tumors with wild 
type p53 it is common for p53 expression to be misregulated. For 
example, proteins that have a part in downregulating p53, such 
as MDM2 and MDMX, are commonly overexpressed in human 
tumors (Momand et al., 1998; Danovi et al., 2004). Furthermore, it 
has been demonstrated that restoration of p53 function can cause 
tumors to regress in vivo (Ventura et al., 2007). The importance of 
p53-signaling in cancer progression, and its therapeutic implica- 
tions, have been investigated in previous mathematical models 
(Gammack et al., 2001; Perfahl et al., 2011), which highlights 
further our study. 

Note that simply removing the limitations on a cell imposed 
by p53 is not enough for it to become cancerous; for a cell to 
become cancerous it must accumulate multiple hallmarks includ- 
ing: self-sufficiency in growth signals, insensitivity to an ti- growth 
signals, limitless replicative potential, sustained angiogenesis, and 
the ability to migrate to other tissues (Hanahan and Weinberg, 
2011). When such traits accumulate in a cell lacking functional 
p53, the probability of a cell becoming cancerous rises (Alberts 
etal, 1994). 

MDM2 

The protein MDM2 is a key player in the regulation of p53 (Bond 
et al., 2005) and it has been found that MDM2 is commonly ampli- 
fied in human cancers (Momand et al., 1998). MDM2 has been 
shown to be an E3 ubiquitin ligase for p53 (Honda et al, 1997). 
This means that MDM2 can mark p53 for degradation by the 
proteasome. As such, amplification of MDM2 leads to reduced 



DNA repair 
and return to cell cycle 



Apoptosis 



I 



GO 



p53 activation — 1 



Potential 
carcinogenesis 



Failure to hi 
the cell cycl 



iiit] 
:le | 



Benotoxic stress 



p53 activation 
failure 



G1 




S 




G2 




M 








t 











FIGURE 1 | Diagram of p53 and the cell cycle, showing possible 
outcomes of stress and p53 activation. 




FIGURE 2 | Relationships between MDMX, MDM2, and p53. MDM2 
inhibits p53 and is promoted by it. MDM2 inhibits itself and this effect is 
reduced by MDMX. MDMX inhibits p53 directly, and is itself inhibited by 
MDM2. 



p53 levels (Haupt et al, 1997; Kubbutat et al, 1997). MDM2 
production is also induced by p53, forming a feedback loop (Barak 
et al, 1993). Figure 2 illustrates the interactions of MDM2 with 
p53. Additionally, MDM2 helps to regulate itself by autoubiquiti- 
nation, meaning it marks itself for degradation by the proteasome 
(Fang et al, 2000). MDM2 possesses a nuclear localization signal, 
which is a structure on the protein that induces the cell to import 
the protein into the cell nucleus (Chen et al., 1995). MDM2 also 
has a cryptic nucleolar localization signal, which flags the protein 
for localization to the nucleolus, but only when MDM2 binding 
to another molecule changes the conformation of the signaling 
region (Lohrum et al., 2000). 

In 2004 several small molecule inhibitors for the p53-MDM2 
interaction were discovered (Vassilev et al., 2004). One of these 
inhibitors, Nutlin-3, was in Phase I clinical trials for retinoblas- 
toma (Secchiero et al., 2011). Nutlins may also have some p53- 
independent effects, and these may be related to MDM2. It has 
been shown in some cell lines that MDM2 is upregulated by 
hypoxia independently of p53 (Gillespie, 2007). Furthermore, it 
has been shown that Nutlin-3 can radio-sensitize hypoxic cells 
that are p53 null, although it has a greater effect on cells with 
wild type p53 (Supiot et al, 2008). Additionally, Nutlin-3 has 
been shown to bind to several anti-apoptotic proteins other than 
MDM2, further complicating any analysis of its effects (Ha et al., 
2011). MDM2 inhibitors bind to the protean competitively and 
occlude the binding site with p53 (Barakat et al., 2010). To the 
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best of our knowledge Nutlins do not alter the autoubiquitination 
properties of MDM2. 

MDMX 

Another important regulator of p53 is MDMX, a homolog of 
MDM2 (Shvarts et al., 1996; Finch et al., 2002). MDMX is com- 
monly overexpressed in tumors, and its upregulation has been 
shown to promote tumor formation (Danovi et al., 2004). Unlike 
MDM2, however, MDMX expression is not induced by DNA dam- 
age (Shvarts et al, 1996). MDMX binds to both MDM2 (Sharp 
et al, 1999) and p53 (Shvarts et al, 1996). MDMX binding 
to MDM2 inhibits MDM2 autoubiquitination (Okamoto et al, 
2009). Furthermore, MDM2 ubiquitinates MDMX (De Graaf 
et al, 2003). The interaction of MDMX and p53 has been shown 
to inhibit p53 activity (Marine et al., 2007). Figure 2 schemati- 
cally depicts the relationships between p53, MDM2, and MDMX. 
MDMX possesses a cryptic nuclear localization signal (LeBron 
et al., 2006), so it can only reach the nucleus while bound to other 
proteins. MDMX is normally located primarily in the cytoplasm 
(Gu et al, 2002). 

Small molecule inhibitors of MDMX have only recently been 
discovered (Reed et al, 2010). Although initial results show some 
efficacy against cancers with upregulated MDMX in cell culture 
(Wang et al., 2011), more work will need to be done to show 
whether or not they will be active in vivo, as well as whether or not 
it is the MDMX interaction or some off-target interaction that is 
causing the effect. 

UPSTREAM REGULATORS 

There are many feedback loops known to affect p53, and the 
behavior of the p53 system is mediated by a number of upstream 
regulators (Harris and Levine, 2005). For example, the protein 
ATM is activated in response to ionizing radiation (Bakkenist and 
Kastan, 2003). Active ATM phosphorylates p53 (Banin etal, 1998), 
MDM2 (Maya et al, 2001), and Chk2 (Matsuoka et al, 2000). A 
related protein, ATR, phosphorylates p53 in response to single 
strand breaks in DNA (Tibbetts et al, 1999). Chk2 along with 
Chkl also phosphorylate p53 (Shieh et al., 2000). These phospho- 
rylations disrupt the ability of MDM2 to affect p53 (Zhang et al., 
1998; Chehab et al, 2000; Maya et al, 2001). 

OTHER FEEDBACKS 

Aside from the MDM2 loop, there are other feedbacks affecting 
p53, although many of these involve also MDM2. The ARF protein 
is known to bind to MDM2 and promote its degradation (Zhang 
et al, 1998). ARF causes both MDM2 and MDMX to be localized 
to the nucleolus (Weber et al., 1999; Jackson et al, 2001). ARF is 
negatively regulated by p53 in a complex manner, thus forming a 
feedback loop (Stott et al, 1998; Lowe and Sherr, 2003). MDM2 
activity becomes enhanced by a feedback in which p53 upregulates 
cyclin G, which then forms a complex with PP2A phosphatase. 
This complex then dephosphorylates MDM2, removing the inhi- 
bition caused by the phosphorylation effect (Harris and Levine, 
2005). The Wipl protein is induced by p53 and is able to mod- 
ify ATM and Chk2, deactivating these proteins, and thus resulting 
in a stronger interaction between p53 and MDM2 (Fiscella et al, 
1997; Fujimoto et al, 2006; Shreeram et al, 2006). Pirh2 has a 



more direct feedback with p53. Like MDM2, Pirh2 and COP1 
both ubiquitinate p53 and are upregulated by p53 (Leng et al, 
2003;Dornan et al, 2004). 

PROTEIN LEVEL OSCILLATIONS? 

Lahavetal. (2004), Geva-Zatorsky etal. (2006), and Geva-Zatorsky 
et al. (2010) all directly observed sustained oscillations of p53 and 
MDM2 levels in the nuclei of individual cells. It is worth noting, 
however, that these single cell studies used MCF-7 cells. MCF-7 
cells were initially used to study p53 because they exhibit wild type 
p53 (Lahav et al, 2004). Unfortunately, the MCF-7 cell line has a 
mutation in an MDM2 intron causing upregulation of MDM2 (Hu 
et al, 2007), lacks ARF (Stott et al, 1998), and possesses amplified 
MDMX (Danovi et al., 2004). Because of this, any assumption that 
any wild type cell would behave similarly to an MCF-7 cell with 
respect to p53 regulation is questionable at best. Unfortunately, 
there are no similar single cell studies of non-tumorigenic cell 
lines at the time of writing this paper. Also of note is the finding 
by Batchelor et al. (2011) that MCF-7 cells respond differently to 
damage induced by ultraviolet light than they do to double-strand 
breaks induced by gamma radiation or radiomimetic drugs. Geva- 
Zatorsky et al. (2006) also pointed out that undamped oscillations 
of p53 levels may appear damped in studies of cell populations 
due to the individual cells falling out of sync with each other. 
Damped oscillations have been observed in populations of non- 
tumorigenic cell lines, for example in entire mice (Hamstra et al, 
2006). 

PREVIOUS MODELING WORK 

A number of models of p53 response to DNA damage have been 
proposed in the past. These models are based on a variety of 
approaches and serve a number of functions. Some basic models 
use built-in time delays on p53 induction of MDM2 transcrip- 
tion, such as some of the models developed by Geva-Zatorsky 
et al. (2006). In contrast, the model presented by Lev Bar-Or et al. 
(2000) used coupled differential equations to create time delay 
effects. There are advantages and disadvantages to each of these 
approaches. In a real cell, proteins are not produced instantly 
in response to a promoter. Both transcription and translation 
processes take time, and transport of the mRNA and the protein to 
the cytoplasm does not happen instantaneously. An explicit time 
delay deals with this problem directly, but may be more difficult 
to analyze than coupled equations. It also adds to the complex- 
ity of any computer algorithm made for stochastic simulations. 
A set of coupled equations, on the other hand, will start to show 
effects of induced protein production in the protein levels instan- 
taneously, but the effect will be very small until some time has 
passed. In a stochastic system the protein levels are quantized and 
instead of instantaneous effects there is simply a small but non- 
zero possibility of instantaneous effects. In both the stochastic and 
deterministic cases adding more steps in the form of more cou- 
pled equations makes the system both more realistic and more 
computationally intensive. Another factor to consider is that p53 
induces the transcription of MDM2 mRNA, and that mRNA is 
active for a time. Because of this, the actual rate of MDM2 pro- 
duction is dependent on a weighted average of past p53 levels 
rather than p53 levels at some specific time in the past. Using a 
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single delayed p53 term to describe MDM2 production is therefore 
problematic. One way around this problem is to use a delay term 
for the production of the MDM2 mRNA rather than the MDM2 
protein, as was done by Cai and Yuan (2009). Ma et al. (2005) 
investigated the number of p53 pulses that occur in response 
to DNA double-strand brakes using a model made from three 
linked modules, simulating DNA repair, ATM activation, and the 
p53-MDM2 feedback loop. Linking together multiple systems like 
this, in particular linking to systems that can be easily perturbed 
experimentally, may be a good way to develop models that are 
straight- forward to test. Batchelor et al. (2008) proposed a model 
based on abstracted signal and inhibitor systems interacting with 
MDM2 as well as active and inactive p53. This model was cre- 
ated to investigate the possible effects of ATM, CHK2, and WIP1 
on p53 behavior. They included an equation for an input signal 
that converted p53 from an inactive form to an active form, and a 
p53 induced inhibitor that reduced the effects of the signal. There 
have also been past efforts to look at stochastic models of the p53 
regulatory system. Cai and Yuan (2009) modeled p53-MDM2 and 
MDMX interactions and analyzed some of the effects of intrin- 
sic noise. Their model has MDM2 mRNA being produced with 
a time delay. It also includes ubiquitinated states of proteins and 
a deubiquitination term, rather than just assuming all ubiquiti- 
nated proteins are degraded. Puszynski et al. (2008) developed a 
complex stochastic model of p53 behavior aimed at showing how 



stochastic effects lead to variability of cell fate in a bistable model. 
Their model includes a cytoplasmic compartment and a nuclear 
compartment, although p53 is not included in their cytoplasmic 
compartment. In addition to the negative feedback of MDM2 and 
p53 they include a positive feedback involving a series of events 
that lead to MDM2 being sequestered in the cytoplasm where it 
can no longer degrade p53. 

Table 1 summarizes the key differences between the models. 
Ultimately, the differences in the models have as much, if not more, 
to do with differences in what the researchers were trying to inves- 
tigate, rather than with differing assumptions about p53 behavior. 

MATERIALS AND METHODS 
THE MODEL 

Since it has been observed that stochastic effects can cause a pop- 
ulation of cells that undergo undamped oscillations to appear as 
if they were undergoing damped oscillations (Lahav et al., 2004; 
Geva-Zatorsky et al, 2006), it is interesting to compare a stochas- 
tic model of cell behavior to a deterministic one. By using both 
stochastic and deterministic versions of the same model it will be 
possible to look at the process of desynchronization between cells, 
which causes oscillations to appear damped, and to search for any 
other effects by which stochasticity could influence the system. As 
we shall see later, further investigation revealed several unexpected 
ways in which stochasticity influenced the system. 



Table 1 | Key features of various models of p53 behavior. 



Model 


Stochasticity 


MDMX 


Compartments 


Time delayed 
equations 


Stress 
signal 


Other notes 


Geva-Zatorsky 
etal. (2006) 












These models do not have 
saturable MDM2 production 


Model 1 


Limited noise 


No 


No 


No 


No 


Linear Model 


Model 2 


Limited noise 


No 


No 


No 


No 




Model 3 


Limited noise 


No 


No 


Yes 


No 


Linear Model 


Model 4 


Limited noise 


No 


No 


No 


No 




Model 5 


Limited noise 


No 


No 


No 


No 


Linear Model 


Model 6 


Limited noise 


No 


No 


Yes 


Yes 




Lev Bar-Or 
etal. (2000) 


None 


No 


No 


No 


Yes 


Stress is abstract and gets 
repaired 


Ma et al. 
(2005) 


In the stress and repair 
modules only 


No 


No 


Yes 


Yes 


Complex stress and repair 
modules 


Batchelor et al. 
(2008) 


No 


No 


No 


Yes 


Yes 


p53 promotes an inhibitor of 
the stress signal 


Cai and Yuan 
(2009) 


Yes 


Yes 


No 


Yes 


No 


Includes phosphorylated 
proteins 


Puszynski 
etal. (2008) 


Yes 


No 


Yes, but not for p53 


No 


Yes 


Includes many other proteins 


Our model 


Stochastic and non-stochastic 
versions were implemented 


No 


Only for MDM2 


No 


No 


Details in Section "Materials 
and Methods" 



Frontiers in Oncology | Molecular and Cellular Oncology 



April 2013 | Volume 3 | Article 64 | 4 



Leenders andTuszynski 



Models of cellular p53 regulation 



In this model p53 induces the transcription of MDM2 mRNA in 
the nucleus; there are three steps between induced transcription 
of MDM2 by p53 and the arrival of MDM2 proteins in the cell 
nucleus. Induced transcription is assumed to be proportional to 
[p53]/(Kp- 8 + [p53] 1,8 ), as was seen in the binding properties 
found by Weinberg et al. (2005). MDM2 mRNA is also produced 
at a basal rate. After being produced in the nucleus, the MDM2 
mRNA proceeds to the cytoplasm, where it is translated and even- 
tually decays. Even though mRNA from MDM2's different pro- 
moter regions are translated at different rates, they are treated as 
one species. Because the two types of mRNA are assumed to decay 
at the same rate, this amounts to absorbing the difference in trans- 
lation rates into the mRNA production rates. Cytoplasmic MDM2 
moves to the nucleus at a constant rate, and all other behaviors that 
cytoplasmic MDM2 could exhibit are ignored in this model. ARF 
was given constant production and degradation rates. Once in the 
nucleus, MDM2 can become bound to ARF, which removes both 
proteins from the system. Additionally, MDM2 autoubiquitinates, 
which is a process that also removes it from the system. Figure 3 
provides a schematic diagram of this system. 

Using the principle of mass-action and the saturable tran- 
scription kinetics mentioned above, the system's behavior can 
be mathematically described in terms of a system of differential 
equations. In addition to all the chemical reactions in Figure 3 
the system of differential equations includes the production and 
degradation of p53, basal transcription of MDM2 mRNA, decay 
of cytoplasmic RNA, decay of ARF, and production of ARF. The 
equations are as follows: 



d [p53] 



dt 

d [RNA nud ear 



kp - h [p53] [MDM2 nuclear ] - d v [p53] 



[p53] L8 

= k m +^2 t1R 18 - fa [RNA nud ear] 



dt **£ 8 + [p53] 



d [RNA cytop i asm i c ] 



dt 



= k 0 [RNA nuc iear] ~ d TC [RNA cytop i asmic ] 



d [MDM2 cytoplasmic ] 



dt 



d [MDM2 nuclea r] 

dt 



d [ARF] 
dt 



= kj [RNA cytop iasmic] 

- k [MDM2 cytop iasmi c ] 
k [MDM2 cytop iasmi c ] 

- d mn IMDM2 2 , 1 

11111 |_ cytoplasmic J 

- h [MDM2 nuclear ] [ARF] 

= h - 4 [ARF] - k 3 [MDM2 nuclear ] [ARF] 



with k v being the production rate of p53, k\ being the rate 
at which MDM2 ubiquitinates p53, and dp being the rate of 
MDM2 -independent p53 degradation. Here, k m is the rate of p53- 
independent MDM2 mRNA production, kj is the maximum rate 
of p53 -dependent MDM2 mRNA production, Kd is the dissocia- 
tion constant for p53 on the MDM2 promoter region, and ko is the 
rate of MDM2 mRNA transport to the nucleus. In the equations 
above, d YC is the decay rate of MDM2 mRNA in the cytoplasm, 
kj is the translation rate for MDM2 mRNA, and k x is the rate of 



Cytoplasm 



1.66*10^/3 

MDM2 • MDM2 RNA 
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FIGURE 3 | A schematic of the model of p53 including MDM2 
sequestration by ARF. The blue boxes denote molecular species in the 
cytoplasm. The yellow boxes indicate molecular species in the nucleus. 
Arrows denote movement between compartments, barred lines indicate 
degradation, and circles indicate inducing production. 



nuclear localization for MDM2. MDM2 autoubiquitination hap- 
pens at the rate d mn and MDM2 binds to ARF at the rate k$. 
Lastly, ARF is produced at the rate k a and degraded at the rate d a . 
The binding properties of p53 and the MDM2 promoter have been 
investigated experimentally by Weinberg et al. (2005), who showed 
that the appropriate Hill coefficient for the Hill function is 1.8. 

A list of the values used for these parameters can be found in 
Table 2. The initial conditions were chosen by letting the system 
run until it settled into a stable limit cycle and then by using 
the values for the time when nuclear MDM2 levels were at a 
maximum. 

Experimental observations of the p53-MDM2 feedback loop 
have found periods of oscillations between 4 and 7h (Geva- 
Zatorsky et al., 2006, 2010). Due to scarcity of experimentally 
verified data, most of parameters in the model were chosen by 
hand in order to produce oscillations with a similar period. Some 
of the parameters were constrained by experimental data. was 
found to be 12.3 nM by Weinberg et al. (2005). Some exper- 
imental results suggested that the half-life for MDM2 mRNA 
should be in the range of 1-2 h (Hsing et al., 2000; Mendrysa 
et al., 2001), so this constrained our choice of the decay rate. 
The MDM2 translation rate, kj, was assumed to be one protein 
per mRNA molecule per minute, approximately the value esti- 
mated by Cai and Yuan (2009). The transport rate for MDM2 
mRNA was constrained to be in the range of 5-40 min, based 
on Mor et al. (2010). The half-life of the ARF protein, d a , was 
chosen to be 6h based on Kuo et al. (2004). Complex for- 
mation rates were assumed to be 6 x 10 _4 /nMs, a reasonable 
rate for protein-protein interactions (Northrup and Erickson, 
1992). It was further assumed that the p53-MDM2 interaction 
would always result in p53 degradation. MDM2 -independent p53 
turnover was assumed to give a half-life of 10 h for the p53 pro- 
tein; this is essentially negligible in this model, but this term was 
included in the model so that a bifurcation value could be cal- 
culated for it. Cytoplasmic volume was assumed to be 1000 |xm 3 
with a nuclear volume of 100 |xm 3 . The values for p53 production, 
ARF production, basal MDM2 mRNA production, p53 induced 
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Table 2 I Parameters used in the model. 



Parameter 


Description 


Value 


Value (alternate expression) 




p53 Production 


0.5 proteins/s 


Q QH 1H— 3/nN/l c 

o.ou x iu /nivi s 




iviuiviz ciepencieni poo GegidGdiion 


y.yuo x iu /s 


r n/ in - 4/ n i\ /| c 
ox iu /nivi s 


rl 

dp 


p53 Decay 


i.yzo x iu /s 


-in h half lifo 

i u n ridiT-iiTe 


u 

Km 


POO 1 1 lUUpUl 1UUI 1 L IVIUIVIZ piUUUCLIUIl 


i c Y in - 3 rma/q 

I.O X IU nlMM/b 


1 RMA nor fififi c 
i ni\i/-\ pt?i OOO b 




poo-uepencieni iviuiviz prociuciion 


I.O X IU /s 


l\ A ovim i i ry-\ r\-f 1 DMA v~\r\r o 

ividximum ot i ni\i/-\ per oob 


AT'r-s 


UlbbUOId UUI 1 CUNbLdNL 


/^U pi ULcll lb 


i z .o n ivi 


Isrs 

KO 


RNA transport from nucleus to cytoplasm 


o.u x iu /b 


14.4 min for half the proteins to move 


rl 


iviuiviz mniMM aecay in cytoplasm 


1 AAA i n-4 /c 
l.^^-^ X IU /s 


1 h OH min half li-fa 

i n zu mm nau-iue 




Transcription rate 


1.66 x 10~ 2 proteins/s 


One protein per RNA per min 


k\ 


Protein transport from cytoplasm to nucleus 


9.0 x 10" 4 /s 


12.4 min for half the proteins to move 




MDM2 autoubiquitination 


1.66 x 10" 7 /s 


2.76 x 10" 9 /nMs 




ARF production 


0.5 proteins/s 


8.30 x 10" 3 /nMs 


da 


ARF decay 


3.209 x 10" 5 /s 


6 h half-life 


k 3 


MDM2-ARF complex formation rate 


9.963 x 10" 6 /s 


6x 10" 4 /nMs 
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FIGURE 4 | p53 and MDM2 oscillating in the deterministic model. p53 is in black, MDM2 is in red 



MDM2 mRNA production, MDM2 nuclear import, and MDM2 
autoubiquitination were unknown. These unknown parameters 
were chosen manually in order to produce oscillations similar to 
the ones observed in experiments on single cells. Although only 
one set of parameters was produced for this model, the choice 
is certainly not unique given the somewhat loose selection crite- 
ria. The model produces oscillations with a period of 6.4 h as can 



be seen in Figure 4. Bifurcation points for the model are listed 
in Table 3. The bifurcation points were found numerically using 
Matlab (MathWorks, Inc.). 

STOCHASTIC SIMULATION ALGORITHM 

The Gillespie algorithm is one of the most commonly used 
methods of stochastic simulation (Gillespie, 1977). The Gillespie 
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Table 3 | Bifurcation points in the deterministic model. 



Parameter Bifurcation value Oscillatory behavior 





0.21 5/s 




Undamped: 0.215 < k p < 1.462 




1.462/s 




Damped: k p < 0.215U k p > 1.462 




2.903 x 10" 


" 6 /s 


Undamped: 2.903 x 10" 6 <k<\< 
1.834 x 10" 5 




1.834 x 10" 


" 5 /s 


Damped: k<\ < 2.903 x 10" 6 U k<\ > 
1.834 x 10" 5 


d P 


4.237 x 10" 


" 4 /s 


Undamped: d p < 4.237 x 10" 4 


km 


2.788 x 10" 


" 3 /s 


Undamped: k m < 2.788 x 10" 3 


k 2 


7.501 x 10" 


3 /s 


Undamped: 7.501 x 10" 3 < k 2 < 0.118 


k 2 


0.118/s 




Damped: k 2 < 7.501 x 10" 3 U k 2 > 0.118 


Kd 


253.083 




Undamped: 253.083 < K D < 1723.058 


Kd 


1723.058 




Damped: K D < 253.083 U K D > 1723.058 


ko 


7.010 x 10~ 


6 /s 


Undamped: 7.010 x 10" 6 <k 0 < 
6.160 x 10" 3 


ko 


6.160 x 10" 


" 3 /S 


Damped: k 0 < 7.010 x 10" 6 U k 0 > 
6.160 x 10" 3 




8.714 x 10" 


" 5 /S 


Undamped: 8.714 x 10" 5 < c/ rc < 
2.704 x 10" 4 


drc 


2.704 x 10" 


" 4 /s 


Damped: c/ rc < 8.714 x 10" 5 U d rc > 
2.704 x 10" 4 


kj 


8.760 x 10" 


- 3 /s 


Undamped: 8.760 x 10" 3 < kj < 
2.936 x 10" 2 


kj 


2.936 x 10" 


" 2 /s 


Damped: kj < 8.760 x 10" 3 U kj > 
2.936 x 10" 2 


k\ 


6.845 x 10" 6 /s 


Undamped: 6.845 x 10" 6 <k\< 
1.559 x 10" 2 


k\ 


1.559 x 10" 


" 2 /s 


Damped: k\ < 6.845 x 10" 6 U k\ > 

1 ceo 1 n— 2 

I. boy X ID 


dmn 


1.251 x 10" 


" 6 /S 


Undamped: c/ mn < 1.251 x 10" 6 


k a 


0.324/s 




Undamped: 0.324 <k a < 0.963 


k a 


0.963/s 




Damped: k a < 0.324 U k a > 0.963 


da 


2.088 x 10" 


" 3 /s 


Undamped: d a < 2.088 x 10" 3 


ks 


5.866 x 10" 


" 6 /S 


Undamped: k 3 > 5.866 x 10" 6 



algorithm has the advantage of being exact, unfortunately, it is also 
computationally expensive. In order to conduct our investigation 
we chose to instead use an approximate simulation, because the 
Gillespie algorithm is too slow for the required complexity and 
number of simulation runs. 

The algorithm we created was based on the concepts of a finite 
difference integrator. In a finite difference integrator a system of 
differential equations is evaluated by first calculating each of the 



derivatives at a point in time, then multiplying them by the time 
step size, and finally updating each of the variables by the corre- 
sponding amount. In our algorithm, rather than being evaluated 
as a single set of derivatives each chemical reaction is evaluated 
separately. When the simulation evaluates a chemical reaction, the 
first step is to use the law of mass -action and the average of the 
current chemical concentrations, and their concentrations after 
the last time the reaction was evaluated, to find an expectation 
value for the number of times the reaction will occur during this 
time step. Next, the expectation value for the number of times 
the reaction will occur is set as the expectation value for a Pois- 
son random number generator and the result is the number of 
times the reaction will actually occur during that time step. This 
gives the algorithm a strong resemblance to the well known tau 
leap method (Gillespie, 2007), in which Poisson random num- 
bers are used in combination with the Gillespie algorithm to 
improve efficacy. In order to improve efficiency while preserv- 
ing accuracy in our algorithm, an adaptive time step is used. The 
program evaluates each reaction 0.5^ times per simulated second, 
with N chosen such that the expectation value for a particular 
evaluation of a reaction is lower than a preset threshold mul- 
tiplied by the quantity of the chemical molecules involved. In 
this way parts of the system that are changing rapidly are eval- 
uated with a low enough time step to prevent numerical errors, 
without needing to waste additional computations on the slower 
reactions. 

Figure 5 shows some examples of individual simulation runs 
for this model. The stochastic nature of the simulation leads to 
a number of interesting differences arising from the desynchro- 
nization of the individual model runs as well as from applying a 
distribution of p53 values into the non-linear function for MDM2 
production. 

RESULTS 

DESYNCHR0NIZATI0N IN GENERAL 

In order to understand how the individual stochastic realizations 
of our model fall out of synch with each other let us first consider 
how stochastic systems may fall out of synchronization in general. 
An experiment averaging protein levels across many cells is analo- 
gous to looking at the average of many runs of a stochastic system. 
As such, it is interesting to consider how aggregate average behav- 
ior differs from the behavior of individual model runs. A given 
run of the stochastic model will not necessarily just be equal to the 
deterministic model plus noise. At any given step the stochastic 
model's variables depend on the values of the variables at the pre- 
vious time step. For a periodic model this will result not only in 
noise moving variables up and down but also in random stepping 
forwards and backwards of the model's phase. Consequently, an 
ensemble of model runs will fall out of synchronization over time. 
Imagine for simplicity a stochastic model based on a deterministic 
model with a variable given by A sin (cot + cp). In the stochastic 
model random chance continuously moves each run in the ensem- 
ble toward or away from the next peak. Considering the central 
limit theorem applied over a large number of runs, one would then 
expect the distribution of timing of the peak in individual runs to 
approach a normal distribution. If all the runs are initialized from 
the same starting point, then the amplitude of the mean will not 
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FIGURE 5 | Examples of time courses in the stochastic model. p53 is in black and MDM2 is in red. 



be A sin (cot + cp) but rather it will be 
1 



J— c 



e 2a 2 sin (oof + cp + (tit')dt' 



will increase proportionally to the square root of time, the stan- 
dard deviation a can be expanded as a^/t> where a is a parameter 
related to the rate of desynchronization. This integral then works 
out to be 



because the timing of each run will be shifted with a Gaussian 
weighting given to the shift. Since the width of the distribution 



Ae~ 



2 sin (oof + cp) 
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FIGURE 6 | Comparison of stochastic and deterministic models. 

(A) Shows the comparison for MDM2 with MDM2 from the 
deterministic model in red and from the mean of 5,000 runs of the 



stochastic model in blue. (B) Shows the comparison for p53 with the 
deterministic model in black and from the mean of 5,000 runs of the 
stochastic model in green. 




FIGURE 7 | (A) Comparison of p53 levels in the deterministic model in black to a curve fitted to it from the function f (f) = a 0 + a : sin (oof + <p0 + 

a 2 sin (2 oof + cp 2 ) in red. (B) Comparison of p53 levels in the stochastic model in black to a curve fitted to it from the function f (f) = a 0 + e~~ a^ sin (oof + <p0 - 

e -4 ^ a 2 sin (2oof + <p 2 ) in red. 



Applying the result above we find that the function will be 



Consider a 2tt periodic function that is integrable on the inter- 
val from —7i to 7t. This function could be expressed as a Fourier changed by desynchronization to become 
series such that 



/ (t) = — + ^ [dn cos (nt) + b n sin (nt)] 



n=l 



or equivalently 



/<->-?+ 



oo 

^a n sin (^nt + + b n sin (nt) j 



2 2 

f (t) = y + ^ |^a„ sin ^nt + + fc M sin (nt)j^ _ZL ^ 



Since the decay is proportional to the square of the frequency, 
any function will rapidly take on the appearance of a single 
decaying sine-function curve as time progresses. 
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DESYNCHRONIZATION IN THE STOCHASTIC MODEL 

The damping caused by desynchronization in the stochastic model 
can be seen in Figure 6. The deterministic and stochastic systems 
can be compared by fitting a curve to the time series for p53. 
Specifically: 

f (t) = oq + a\ sin (oof + cpi) + sin (2oof + 92) 
for the deterministic model, and 

f (t) = oq + e 2 fli sin (oof + cpi) + e 2 fl2 s in (2oof + 92) 

for the stochastic model. Table 4 lists the parameter estimates for 
the deterministic model as well as 95% confidence intervals for the 
stochastic model. Figure 7 shows graphs of the functions and their 
best fits. The best fit was determined by using least squares regres- 
sion on the mean p53 values from 5,000 instances of the stochastic 



Table 4 | Comparisons of the parameters found when fitting the 
deterministic model's p53 levels to the function / (t) = ao + 

a\ sin (of + cpi) + a 2 sin (2wf + (p 2 ) and the stochastic model's p53 

ot 2 t 

levels to the function / (t) = ao + e~ ~ a\ sin (oof + (pi) + 

_ 4ot 2 t _ 

e 2 #2 sin (2(ot + (p2). 





Parameters 


Parameters 


Lower bound 


Upper bound 




fitted to 


fitted to 


for stochastic 


for stochastic 




deterministic 


stochastic 


parameters 


parameters 




model 


model 






a 


N/A 


21.8/s 1/2 


21.2/s 1/2 


22.5/s 1 / 2 


00 


2.73 x 10" 4 /s 


2.63 x 10" 4 /s 


2.62 x 10- 4 /s 


2.64 x 10" 4 /s 


A 0 


332 
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A, 


-348 


-396 


-406 
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1.21 


1.40 


1.38 
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A 2 


105 


136 


-136 


144 
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-1.16 


-36.6 


13.6 
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Time (h) 

FIGURE 8 | Comparison of nuclear MDM2 levels in the stochastic 
model in black to a curve fitted to it from the function 

a 2 t 4a 2 t 

f (t) = ao + e~ t flj sin (wf + (p x ) + e~ ~ a 2 sin (2(*)f + (p 2 ) in red. 



model. The upper and lower bounds were found by using boot- 
strapping on the 5,000 instances that were used to compute the 
best fit. The 95% confidence intervals for the amplitude and phase 
of the second sine curve ended up being very large due to the curve 
fitting function jumping between local minima. To ensure that the 
algorithm was being run at a high enough numerical precision, 
an additional 5,000 instances were generated with the acceptable 
error parameter in the code selected to equal 10 times the value 
used in this analysis. The resulting new confidence intervals were 
compared to the ones from the higher accuracy runs. In all cases 
significant overlap of the intervals was found, suggesting that the 
acceptable error was set low enough in the high accuracy runs to 
result in only negligible deviations from an exact solution. 

The differences between the stochastic model's behavior and the 
deterministic model's behavior are statistically significant. Most 
striking is that the frequency of the oscillations was changed by 
stochastic effects. The same analysis has been done on nuclear 



Table 5 | Comparisons of the parameters found when fitting the 
deterministic model's nuclear MDM2 levels to the function / (t) = a 0 + 

a\ sin (cat + (pi) + a 2 sin (2wf + (p 2 ) and the stochastic model's nuclear 
MDM2 levels to the function / (t) = oq + e~~ a\ sin (of + (pi) + 
e 2 a2 sin (2wt + 92). 
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FIGURE 9 | A comparison of the function between the 

function applied to mean p53 values in black and the mean of the 
function when applied to the distribution of p53 values in red. 
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MDM2 levels, which can be seen in Figure 8 and Table 5. The 
discrepancy between the fitted curve for MDM2 levels and the 
levels from the simulation hints at another difference between sto- 
chastic and deterministic systems, which will be discussed below. 
It is also worth noting that this stochastic model only consid- 
ers the differences between cells due to noise in a few chemical 
reactions. In a real cell there would be many more factors con- 
tributing to desynchronization. Even simply adding mRNA for the 
p53 and ARF included in this model raises the desynchronization 
parameter a from 21.8 to 23.5 s _1/2 (a mean of 30 mRNA mol- 
ecules was used for this simulation). Additionally, differences in 



cell volume would increase desynchronization by altering protein 
concentrations between cells. 

CHANGES DUE TO NON-LINEAR EFFECTS 

The mean of a stochastic ensemble for the stochastic model devi- 
ates from the deterministic model not just from desynchronization 
but also due to non-linear effects. For a non-linear function 
applied to a distribution of inputs, the mean of the function 
will not necessarily be equal to the function of the mean. In 
other words, as is well known in statistics, the following is usu- 
ally true: <f(x)> ^f(<x>), unless f is a linear function of x. 
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FIGURE 10 | Comparison of stochastic and deterministic models when 
p53 production is near the lower bifurcation point. (A) Shows the 
comparison for MDM2 with MDM2 from the deterministic model in red and 



from the mean of 5,000 runs of the stochastic model in blue. (B) Shows the 
comparison for p53 with the deterministic model in black and from the mean 
of 5,000 runs of the stochastic model in green. 




40 
Time (h) 



B 8- 



o 
o - 

LO 



o 

u. o 

CD O - 

-Q CO 

E 



o 
o - 

C\J 



o 

o - 




20 



1 

40 
Time (h) 



60 



80 



FIGURE 11 | Comparison of stochastic and deterministic models when 
p53 production is near the upper bifurcation point. (A) Shows the 
comparison for MDM2 with MDM2 from the deterministic model in red and 



from the mean of 5,000 runs of the stochastic model in blue. (B) Shows the 
comparison for p53 with the deterministic model in black and from the mean 
of 5,000 runs of the stochastic model in green. 
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FIGURE 12 | Comparison of stochastic and deterministic models when 
the p53 production value is past the upper bifurcation point. (A) Shows 
the comparison for MDM2 with MDM2 from the deterministic model in red 

Production of MDM2 mRNA in this model is clearly non-linear 
because it is proportional to /(p53) = ^ p53 ^ — j-g-. Figure 9 

fc^ + [p53] ' 

compares the function of the mean to the mean of the func- 
tion for this case. Mean MDM2 values in the stochastic model 
are determined by </[p53)> (the red curve in Figure 9) which 
has a different amplitude then /(<p53>) (the black curve in 
Figure 9). This discrepancy causes the behavior of the system to 
change relative to the deterministic case, which only has mean 
p53 values. This is also the most likely source of the discrep- 
ancy between the fitted curve in Figure 8 and the actual levels 
of MDM2. With production that behaves differently, the initial 
conditions in the simulation would not have represented a point 
on the limit cycle for MDM2 levels. As a consequence, the system 
would have been moving toward the limit cycle at the same time 
as it was desynchronizing. The simple fitted curve cannot possi- 
bly account for this, which is why it did not fit well. p53 levels 
would also have been affected by this but this does not seem to 
have been a large enough effect to be readily noticeable on the 
graph. 

Although the effect on the amplitude of the oscillations with 
the original parameters was relatively small, amounting to approx- 
imately 5%, the non-linear effects can be larger in other situations. 
Consider the case when the p53 production rate is set near to the 
lower bifurcation point, as shown in Figure 10. In this case the 
mean level of MDM2 from the stochastic model ends up being 
higher than the maximum amplitude of the oscillations in the 
deterministic model. A similar phenomenon occurs when p53 
production is near the upper bifurcation point as is shown in 
Figure 11. 

EXCURSIONS FROM THE MEAN 

Stochastic effects continue to play an interesting role in the sys- 
tem's behavior even as we move past the upper bifurcation point, 




Time (h) 

and from the mean of 5,000 runs of the stochastic model in blue. (B) Shows 
the comparison for p53 with the deterministic model in black and from the 
mean of 5,000 runs of the stochastic model in green. 
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FIGURE 13 | Comparison of the function j^j^o between the 
function applied to mean p53 values in black and the mean of the 
function when applied to the distribution of p53 values in red. 



so that the deterministic model exhibits damped oscillations. For 
Figures 12-14, p53 production was set to 1.6, putting the system 
into the realm of damped oscillations. In Figure 12 we can see 
that as the oscillations decay, the MDM2 levels settle in at a value 
significantly higher in the stochastic model than the determinis- 
tic one. From Figure 13 we can see that the non-linear effects of 
variable p53 levels are still altering behavior, but something more 
is occurring this time. In Figure 12B we see that mean p53 levels 
are settling in at a level higher in the stochastic model than in 
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the deterministic one. This seems strange in light of the higher 
MDM2 levels but Figure 14 shows the reason. The stochastic 
nature of the system is sufficient to cause significant excursions 
from the mean even though the oscillations should be decaying. 
Some of the oscillations that occur later on are even larger than 
the initial pulse. Similar behavior has been observed in other sto- 
chastic models such as the one presented in McKane and Newman 
(2005), but has not been previously observed in a model of the 
p53 system. 



DISCUSSION 

The stochastic work we present in this paper differs from previ- 
ous modeling efforts in that its goal is primarily to compare the 
behavior of stochastic and deterministic realizations of the same 
model. This requires only a simple model; therefore much of the 
complexity of the p53 system can be ignored. Since the model 
presented in this work is not aimed at addressing DNA repair, or 
dealing with the problem of variable damage being done, it does 
not include such systems. The model presented here also differs 
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FIGURE 14 | Examples of individual stochastic realizations when p53 production is past the upper bifurcation point. p53 is in black MDM2 is in red. 
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from the previous models in a few other ways. Unlike in other 
models, MDM2 autoubiquitination was assumed to happen at a 
rate proportional to the square of MDM2 concentration. Given 
that MDM2 forms heterodimers with MDMX (Sharp et al., 1999), 
that MDMX inhibits MDM2 autoubiquitination (Okamoto et al, 
2009), and that MDM2 ubiquitinates MDMX (De Graaf et al, 
2003), it seems likely that one MDM2 molecule is ubiquitinating 
a second MDM2 molecule. 

The work on the deterministic and stochastic models presented 
here demonstrates that the effects of stochasticity on the behavior 
of genetic regulatory networks cannot be dismissed without care- 
ful consideration. In our system stochastic effects altered every 
aspect of system behavior. In addition to desynchronization lead- 
ing to the appearance of decaying oscillations, the amount of 
MDM2 in the system increased and the period of the oscilla- 
tions changed. The changes in MDM2 levels became more obvious 
when p53 production was near bifurcation points. When the sys- 
tem was put into a state with decaying oscillations, the quantity 
of MDM2 still remained above that in the deterministic model, 
showing that stochasticity still alters behavior as the system is near 
a steady state. Furthermore, stochastic systems will not necessar- 
ily undergo damped oscillations even when assigned parameters 
that would cause damped oscillations in a deterministic system. 
Instead, they may show sporadic oscillation-like excursions from 
the mean behavior. It would seem then that even for cells in a 
steady state, the distribution of protein levels across a popula- 
tion and over time could wreak havoc with attempts to model cell 
behavior. This has implications for researchers wishing to model 
cell-level processes, as systematic errors could occur in determin- 
istic models with no obvious way to compensate for them. As 
computers and algorithms improve, it may be the case that simply 
moving to stochastic modeling of cell populations will become the 
most practical solution. 

The demonstration that stochasticity can be relevant is very 
general, but it was also shown that the magnitude of the effects 
could vary significantly between systems. The effect on mean pro- 
tein levels could be around 5%, as in the original parameter set, 
or around 50%, as in some of the parameter sets with differing 
p53 levels. The obvious way to experimentally test the relevance 



of stochasticity on any given system is by comparing data from 
cell populations to data from individual cells. Such experimental 
comparisons were, after all, the inspiration for investigating sto- 
chasticity in this system in the first place. The difference between 
a stochastic model and a deterministic one with different parame- 
ters are not likely to be obvious from population data, even if the 
effects of stochasticity are expected to be large. Testing the details of 
stochastic models will require investigating the behavior of indi- 
vidual cells. Of course, stochasticity is not the only factor that 
could drive individual cells to different behaviors. Factors such as 
differences in cell size, different cell cycle stages, and non-uniform 
distributions of components in cell culture medium could all alter 
behavior on the scale of single cells. Untangling these effects is 
potentially a fruitful area for future research. 

A valuable way of expanding the utility of our model would 
be to link it to other models of related processes. The DNA 
repair and damage detection modules in Ma et al. (2005) would 
be a good example of this. Once one system is sufficiently well 
understood, it would be possible to begin analyzing how alter- 
ing it changes connected systems, or conversely, how changing 
connected systems alters it. This could allow one to study down- 
stream drug effects. For that kind of work it would likely be 
best to start as far upstream as possible, in order to facilitate 
the experimental control of inputs. For example, for the p53 sys- 
tem it would make sense to start with a model that quantifies the 
damage ionizing radiation causes to DNA and other cellular sys- 
tems, because the level of radiation a cell is exposed to can be 
controlled in the lab. Then, once that is modeled accurately, one 
could study the DNA damage detection systems, and finally the p53 
response. Repeating this process for other forms of damage, like 
for example ultraviolet light, could bring insight into the system's 
behavior. 
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